Helicity dynamics in stratified turbulence in the absence of forcing 
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A numerical study of decaying stably-stratified flows is performed. Relatively high stratification 
(Froude number ~ fO~^ — fO~^), and moderate Reynolds (-Re) numbers [Re ~ 3 — 6 ■ fO'^) are con- 
sidered, and a particular emphasis is placed on the role of helicity (velocity- vorticity correlations), 
which is not an invariant of the non-dissipative equations. The problem is tackled by integrating the 
Boussinesq equations in a periodic cubical domain using different initial conditions: a non-helical 
Taylor-Green (TG) flow, a fully helical Beltrami (ABC) flow, and random flows with a tunable he- 
licity. We show that for stratified ABC flows helicity undergoes a substantially slower decay than for 
unstratified ABC flows. This fact is likely associated to the combined effect of stratification and large 
scale coherent structures. Indeed, when the latter are missing, as in random flows, helicity is rapidly 
destroyed by the onset of gravitational waves. A type of large-scale dissipative "cyclostrophic" bal- 
ance can be invoked to explain this behavior. No production of helicity is observed, contrary to 
the case of rotating and stratified flows. When helicity survives in the system it strongly affects 
the temporal energy decay and the energy distribution among Fourier modes. We discover in fact 
that the decay rate of energy for stratified helical flows is much slower than for stratified non-helical 
flows and can be considered with a phenomenological model in a way similar to what is done for 
unstratified rotating flows. We also show that helicity, when strong, has a measurable effect on the 
Fourier spectra, in particular at scales larger than the buoyancy scale for which it displays a rather 
flat scaling associated with vertical shear. 

PACS numbers: 47.55. Hd, 47.27.-i, 47.35. Bb, 47.27.ek 



I. INTRODUCTION 

Stratified flows, characterized by the Brunt- Vai'ssala 
frequency, TV, are encountered in many astrophysical 
and geophysical settings, with the ratio of the grav- 
ity wave period Ty/ = '^/N to the eddy turn-over time 
tnl = Lq/Uq being rather small. This ratio is the 
Froude number, Fr = Uo/[NLo], where Uq and Lq are 
the characteristic velocity and length scales of the fluid. 
For example, in the atmosphere, Froude numbers of the 
order of 10~^ or 10"'^ are encountered, whereas in the 
oceans they can be ten times smaller. Such flows, due to 
their deep connection to and their ubiquity in our envi- 
ronment, have been studied and reviewed extensively for 
their scaling and statistical properties [l|, [Ij and for their 
physical structures 0] , as well as in the context of wind 
energy. 

With Fr << 1, the waves are fast compared to turbu- 
lent eddies; they can lead to the formation of stratified 
layers with strong vertical gradients, to intensified mixing 
[j| , in particular in the ocean [sj , and they can dominate 
the dynamics, at least at scales larger than the so-called 
buoyancy scale Lb at which Fr = 1, namely Lt ~ Uq/N. 
This has led to the development of a suite of approxima- 
tions, simpler and easier to integrate than the primitive 
equations, using either asymptotic expansions (see, e .g., 
[ill 0412 ) > two-point closures of turbulence [l3l- 151 . 



well as weak turbulence statistical approaches 16-^. 
As one considers the Froude number at scale £ with U£ 



the velocity at that scale, the Brunt- Vai'ssala frequency 
can be considered as constant vis d vis scale dependency, 
whereas the eddy turn-over time, £/u(, decreases with 
scale. As one moves to smaller scales, non-linear interac- 
tions arc more effective and one can define the scale £ at 
which tml{£) = Tw{£)\ this scale is called the Ozmidov 
scale Loz to be defined below, and one expects that at 
scales smaller than Loz^ isotropy and classical turbulent 
scaling recover. 

Gravity waves (and inertia-gravity waves when rota- 
tion is included) are essential to understand the dynam- 
ics of large scales, the very scales where most of the 
energy resides. Under the assumption of stationarity, 
weak nonlinearitics, weak dissipation and negligible forc- 
ing, one then obtains a balance which involves the pres- 
sure gradient and the gravity (together with the Cori- 
olis force in the presence of rotation). Such a concept 
of balance has proven useful time and again in meteo- 
rology and oceanography, and many variants and issues 
have been considered to study the dynamics of synop- 
tic scales. However, gravity waves couple nonlinearly 
on slow time scales, and undergo wave steepening and 
breaking [2[ , through resonant interactions as already de- 
scribed early on through the weak turbulence formalism 
[igl . [20| . This leads for example to frontogenesis [111, , 
as observed in the upper atmosphere or in the oceans. 
Overturning is also present in direct numerical simula- 
tions (DNS) of stratified turbulence at high resolution 
and high Reynolds number [2^, [2^. Thus, stratified tur- 
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bulence is known to play a role in the vertical mixing of 
deep layers in the oceans, as observed for example in the 
central Pacific ocean as well as in river bends and 
estuaries [29|. 



On the other hand, the role of helicity, the correlation 
between the velocity field u and its curl, the vorticity 
w, has not received much attention outside the realm of 
astrophysics when considering the growth of large-scale 
magnetic fields [27| • It is estimated routinely in the atmo- 
sphere, in conjunction with CAPE (convective available 
potential energy) in order to gauge the possibility of su- 
percell convective storms to become strong [2^, and it 
may be a factor to take into account in the formation of 
hurricanes [1^ (see also [13, HH ) . Taking the curl of the 
geostrophic balance (or in the absence of rotation, the 
so-called cyclostrophic balance) eliminates the pressure 
term, and thus also the thermodynamics that could en- 
ter through pressure gradient coupling with density vari- 
ations and, for example, moisture. Thus, considering the 
dynamical role of helicity in stratified turbulence may al- 
low for a focus on the dynamics of the atmospheric wind 
or oceanic current through decoupling from the pressure 
field. Note that cyclostrophic balance implies curvature 
of trajectories that is not due to rotation (which is ne- 
glected here), as is also the case at low latitudes or at 
small scales, for example in tornadoes [33 |. 



Specifically, in this work wc consider the effect of he- 
licity in freely decaying stably stratified turbulence using 
direct numerical simulations. We show that helicity un- 
dergoes a very slow decay (when compared with the en- 
ergy), and that cyclostrophic balance at large scales can 
play an important role in determining the decay of the 
fiow. 



In Sec. |TT] we recall the incompressible Boussincsq 
equations for a stably stratified fiow and wc list some 
characteristic length scales, dimcnsionlcss number, and 
energies, relevant to the problem under consideration. 
We shortly present the numerical method used to inte- 
grate the equations, and we describe the different ini- 
tial conditions imposed, reporting, in Table HI the initial 
parameters chosen for our calculations. In Sec. IIIII we 
mainly discuss the behavior of helicity for different ini- 
tial conditions, we relate it to the velocity distribution 
in the horizontal and vertical directions and wc inter- 
pret these results proposing a type of large-scale dissipa- 
tive "cyclostrophic" balance. In what follows we study 
how the combined effects of different residual helicity val- 
ues, different initial conditions and different Froude and 
Reynolds numbers affect the temporal decay of the total 
energy and the Fourier spectra for the total and potential 
energy and for the absolute value of helicity. Finally our 
results are summarized in IIVI 



II. EQUATIONS AND METHODOLOGY 
A. Mathematical model and relevant parameters 

Stably stratified turbulence is studied here by means of 
the incompressible Boussinesq equations. In the presence 
of gravity g and in the absence of forcing, with 9 being the 
potential temperature fluctuations in units of velocity, 
the equations read: 



+ u • Vu = -VP - NO 
dtO +u-Ve ^ Nw + kAO , 
V-u = ; 



(1) 
(2) 
(3) 



w is the vertical {z) component of the velocity in the di- 
rection of gravity, P is the pressure, v the viscosity, and k 
the diffusivity; in this paper we take a unit Prandtl num- 
ber, Pr = v/k ^ 1. The square Brunt- Vaisala frequency 
is given by iV^ = ^(.9/^* ){dO /dz), where d9 /dz is the 
imposed background stratification, which is assumed to 
be linear. No filtering on the small scales {e.g., hyper- 
viscosity) is applied, a normal Laplacian being used for 
diffusion. 

In the ideal case (with v = k ~ Q), the Boussinesq 
equations conserve the total (kinetic plus potential) en- 
ergy, 



and the point-wise potential vorticity 

PV = -Nuj^ +LJ-VB , 

which is a material invariant. For strong waves (or weak 
nonlinearities). A'' >> 1, and the linearized version of PV 
leads to a quadratic invariant, (PV^), which is preserved 
by the Fourier truncation of the numerical algorithm and 
which is proportional to the so-called vertical enstrophy. 

The volume integrated helicity, 
H = (u-uj) 

is a topological quantity that characterizes the amount 
of "twist" and linkage in the flow |33|. In the 3D Euler 
equations, in the absence of stratification and of dissipa- 
tion, helicity is an invariant, as is the kinetic energy, but 
in the above equations it is not. We nevertheless retain 
this diagnostic in order to make connections to isotropic 
and homogeneous turbulence, and to rotating fiows for 
which helicity is invariant as well. 

Relative helicity ay measures the degree of alignment 
between velocity and vorticity, and as such, it can be 
viewed as a proxy measure of the amount of non-linearity 
in the flow, since ay ^ ±1 implies a Lamb vector — u • 
Vu + V[u'^/2] = u X cj w 0. It can be defined as 



H 



ay 



\/ Ey Z\ 



(4) 
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where Zy = (w^) is the so-called kinetic enstrophy, pro- 
portional to the kinetic energy dissipation when viscosity 
is restored. One can define similarly a potential enstro- 
phy Zp = (|V6'p), associated with the dissipation of 
potential energy Ep when k 7^ 0. 

Three characteristic length scales can be identified for 
this problem: (z) the Kolmogorov length scale where 
dissipation prevails for a Kolmogorov spectrum, (ii) the 
Ozmidov length scale, Lqz, and (ivi) the buoyancy length 
scale L^,. For the initial conditions, these length scales 
are defined respectively as 



n = 



1/4 



in 



3/2 



u = 



2nUn 



N 



where Uq is the initial rms velocity, Lq — 27r/fco is the 
initial integral scale, fcg being the wave number where 
the initial excitation is centered, and ey = dEy/dt w 
Uq/Lq has been used as an estimate of the kinetic energy 
dissipation rate (under the assumption that the Ozmidov 
scale is resolved). 

Three relevant dimensionless parameters can also be 
identified: (i) the Reynolds number Re = UqLq/v, [ii) 
the above mentioned Froude number, Fr = Uq/ (LqN), 
and {Hi) the buoyancy Reynolds number TZ = ReFr^. 

We finally define the so called reduced spectra E as 
a function of the wavenumbers k, ki_ and fc|| ; they are 
respectively k = |k|, fc^ = |kj^| = |k|sin9, with O the 
co-latitude in Fourier space with respect to the vertical 
axis of unit vector z. while k\\ refers to the component 
of k in the z direction. If [/y(k) is the velocity auto- 
correlation function in Fourier space, we name its trace 
C/(k). Under the assumption of homogeneity we define 
the axisymmetric energy spectrum as 

e(|ki|,fc,|) = y"[/(k)|k|sined0 = e(|k|,e), (5) 

where (j) is the longitude with respect to the x axis; then 
we can define |3J1 



(6) 
(7) 
(8) 



Similar definitions hold for /i(|k|,8), the axisymmetric 
helicity spectrum which is based on the antisymmet- 
ric part of the velocity correlation tensor [s^]- From 
/i(|k|,0) we then derive the definitions for H{k), H{k±) 
and iJ(/c||). Similar definitions hold as well for the po- 
tential energy distribution. 



E{k^) = J e(|kx|,fc||)&||, 
J edkiUii)*^, 
E{k) = / e(|k|,e)|k|de. 



B. Numerical model 

The numerical simulations are carried out using the 
Geophysical High-Order Suite for Turbulence (GHOST) 
code [3^ [3^ . It is a pseudo-spectral framework available 
to the community and hosts a variety of partial differen- 
tial equation (PDF) solvers optimized for studying tur- 
bulence in each of these PDF systems, in two-dimensional 
(2D) and three-dimensional (3D) geometries, for neutral 
and conducting fluids (MHD and Hall-MHD), and with 
solid-body rotation and stratification in the Boussinesq 
approximation. GHOST also includes a passive scalar 
solver, as well as surface quasi-geostrophic and shallow 
water solvers, and several types of sub-grid scale models. 
In the 3D case, the grid is a cubic [0, 27r] '^-periodic box, 
and with a 2"''-order explicit Runge-Kutta (RK) time 
stepping scheme; tests were conducted using a fourth or- 
der explicit RK scheme that showed neither qualitative 
nor quantitative differences when compared with the sec- 
ond order scheme at the resolutions used in this paper. 
De-aliasing is done by using a standard 2/3 rule. The im- 
portance of performing computations for this system in 
a cubic box with equal grid spacing in the horizontal and 
vertical directions was emphasized and demonstrated in 
[33 ■ We also use an explicit time stepping method be- 
cause of the necessity to resolve all time scales. For the 
code, a hybrid MPI/OpenMP parallelization methodol- 
ogy was developed [3^, with tests up to 98000+ compute 
cores on grids of up to 6144^ points; the code also has a 
third level of parallelization with the recent addition of 
GPU/accelerator support for the FFTs. 



C. Initial conditions 

A variety of runs were conducted at moderate Reynolds 
number in a 3D cubic domain, with isotropic discretiza- 
tion of nx = ny ~ nz = n = 256, or n = 512 points. The 
database of the runs is presented in Tabled with param- 
eters characterizing the initial conditions we employed. 
Other runs were performed that confirm the conclusions 
of this work but which, for the sake of brevity, are not 
described in detail. 

Three sets of initial velocities were investigated: 
Taylor-Green (TG), Arnold-Bcltrami-Childress (ABC), 
and random isotropic. The TG velocity field is given 
by: 

wj'' = Vq'~^ sin(A:oa:;) cos(fcoy) cos(A:oz) , 



TG 
..TG 



cos(A;oa;) sm{koy) cos{koz) , 



, 



(9) 



and is globally non-helical (cry = at t = 0). 

The ABC initial condition for the velocity is specified 
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Initial v 


ko 


V 


Re 


Fr 


7^ 


n 




U 


A 


1. TG 


2-3 


2.4e-04 


K. 3000 


0.022 


1.452 


0.0041 


0.0083 


0.35 


0.0246 


2. TGFr/2 


2-3 


2.4e-04 


« 3000 


0.011 


0.363 


0.0041 


0.0074 


0.175 


0.0246 


3. ABC2C 


3-4 


3.0e-04 


« 3000 


0.022 


1.452 


0.0044 


0.0059 


0.25 


0.0246 


4. ABC2C2i?e 


3-4 


1.5e-04 


« 6000 


0.022 


2.904 


0.0026 


0.0059 


0.25 


0.0123 


5. ABC2C2J'r 


3-4 


1.5e-04 


« 6000 


0.044 


11.616 


0.0026 


0.0168 


0.50 


0.0123 


6. ABCAfO 


3-4 


3.0e-04 


K. 3000 


oo 


oo 


0.0044 


oo 


oo 


0.0123 


7. ABC 


3-4 


3.0e-04 


« 3000 


0.022 


1.452 


0.0044 


0.0059 


0.25 


0.0246 


8. ABC2i?e 


3-4 


1.5e-04 


K. 6000 


0.022 


2.904 


0.0026 


0.0059 


0.25 


0.0123 


9. ABC2Fr 


3-4 


1.5e-04 


« 6000 


0.044 


11.616 


0.0026 


0.0168 


0.5 


0.0123 


10. ABC4Fr 


3-4 


3.0e-04 


« 6000 


0.088 


46.464 


0.0026 


0.0474 


1.0 


0.0246 


11. ABCFr-/2 


3-4 


3.0e-04 


^ 3000 


0.011 


0.363 


0.0044 


0.0021 


0.125 


0.0246 


12. RNDiVO 


3-4 


3.0e-04 


^ 3000 


oo 


oo 


0.0044 


oo 


oo 


0.0123 


13. RND 


3-4 


3.0e-04 


7^ 3000 


0.022 


1.452 


0.0044 


0.0059 


0.25 


0.0246 


14. RND2J?e 


3-4 


1.5e-04 


^ 6000 


0.022 


2.904 


0.0026 


0.0059 


0.25 


0.0123 


15. RND2f r 


3-4 


1.5e-04 


6000 


0.044 


11.616 


0.0026 


0.0168 


0.50 


0.0123 


16. RND4f r 


3-4 


1.5e-04 


^ 6000 


0.088 


46.464 


0.0026 


0.0474 


1.0 


0.0123 


17. RNDfc2 


2 


5.24e-04 


3000 


0.022 


1.452 


0.0078 


0.0104 


0.437 


0.0246 



TABLE I. List of the runs, computed on grids with either n=256 or n= 512. The name of the runs ("Initial v") indicates 
the initial velocity field, and also summarizes some important properties (e.g., whether the Reynolds or Froude numbers were 
changed with respect to other runs with similar initial conditions); see the text for more details. All these calculations have 
?7o = 0.5, and dynamical parameters are evaluated at t=0; fco is the wavenumber range of the initial condition, v is the viscosity 
and Re the Reynolds number; Fr is the Froude number, 71. the buoyancy Reynolds number, 77 the Kolmogorov (dissipative) 
length scale, Loz and Lh are the Ozmidov and buoyancy length-scales, and A = 2-k /{n— 1) is the grid resolution. Note that in 
general the Ozmidov length scale is not resolved, except for runs 5, 9 and 14, and that the buoyancy Reynolds number is of 
order unity except for the same three runs. 



as: 



.ABC 

..ABC 
V 

..ABC 



[Bcos{koy) + Csin(fcoz)] , 
[Ccos(fco2^) + ^sin(fcoa;)] , 



..ABC 



[Acos(fcoa-) + B siii(fco2/)] 



(10) 



As a superposition of Beltrami vortices (for which u = 
iw/fco), it is fully helical. Another approximately Bel- 
trami initial condition was taken, namely the v^^'^ writ- 
ten above with v^^'~^ = at t = 0, in order to be able 
to compare with the evolution of the TG flow which also 
has Vz = initially. This second helical initial condition 
is called below ABC2C (2 components). The simulations 
were started using a superposition of these flows (either 
TG, or ABC, or ABC2C) for the wavenumber interval 
reported in Table H] 

We also examined random isotropic initial conditions 
(RND), which use the method described in [sHi to gen- 
erate a helical flow; it does so by correlating the velocity 
and vorticity in such a way that a specific amount of 
relative helicity can be chosen at i = 0. 



In the expressions (O and pUj) . we choose 



„TG 



.,ABC 



1, and, for all the flows, the components of the 
initial velocity are multiplied by a factor /o, such that 
the initial kinetic energy per unit volume is normalized 
to U^: 



V 



{vl+vl + vl)dV 



(11) 



Given Eq. pT|) and the definition of the relative helic- 
ity, Eq. dH), the ABC flow has av{t = 0) = 1, while the 



ABC2C flow displays avit = 0) = /[C^ + {A^+B^)/2l 
which, for the choice A = 0.9, B = 1, C = 1.1, gives 
avit = 0) « 0.57. The choice of A ^ B ^ C allows for 
a more rapid development of the flow and for turbulence 
to set in, by breaking more efficiently the symmetries 
present in the A = B = C case. 

The dimensionless Reynolds number Re, for the TG, 
ABC and random flows, are defined as 



Retg 



Uo2tt 
ko\/3i' 



ReABC.RND 



[/o27r 



The \/3 factor in the definition for the TG flow comes 
from the fact that the individual Fourier shells are pop- 
ulated in different ways, for a given ko, for the TG vs. 
ABC and random flows (e.g., [39|). 

In each run the initial value of the potential temper- 
ature is 9 {t = 0) = 0, so that buoyancy fluctuations 
entirely develop from the initial velocity fleld. Note that 
initial conditions need not be balanced since the time 
step for these runs is small enough to resolve the grav- 
ity waves that develop in each case. A reason for using 
balanced initial conditions is related to the fact that one 
may want to investigate the development of imbalanced 
motions to start as a result of the turbulent mixing that 
develops, rather than by choice of initial conditions. On 
the other hand, the non-linear terms will eventually de- 
velop these unbalanced motions as nonlinear coiling take 
place before and up to the peak of enstrophy. 
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Initial v 


[Ep/Et]* 






cr" 


1. TG 


0.01 


0.0 


0.0 


0.0 


6. ABCAfO 


0.00 


0.99 


0.13 


0.103 


7. ABC 


0.15 


0.99 


0.22 


0.76 


12. RNDiVO 





1.00 


0.13 


0.08 


13. RND 


0.18 


1.00 


0.02 


0.11 


15. RND2fr 


0.18 


1.00 


0.04 


0.13 


17. RNDfc2 


0.17 


1.00 


0.11 


0.70 



TABLE II. Some relevant parameters for a selected set of the runs listed in Table [T] The ratio of potential to kinetic energy is 
given at peak of enstropliy t* , while the relative helicity is reported aX t — Q (cr^), at the peak of the enstrophy (cry), and at 
i = 20 (cry*). Note the increase of relative helicity over time in all stratified cases, after an initial startup phase. 



III. RESULTS 
A. The temporal decay of helicity 

Here we examine the temporal behavior of the flows 
whose initial conditions are given in Table HI For some of 
the simulations, the ratio between the potential energy 
and the total energy at the peak of the enstrophy, the 
value of the relative helicity at the initial time Oy, at 
the peak of the enstrophy cry, and at t = 20, cry*, arc 
reported in Table HIl 

A hint as to the outcome one may expect from the 
time evolution of helicity is obtained by considering the 
dynamical equations. By taking the inner product of Eq. 
(1) with Lj and volume-averaging, it can be shown that 
the time derivative of the total helicity is 



dH 
~dt 



= -2N (9 uj^ - 2vZ 



H, 



Zh = • V X w) 



(12) 



with Zf j th e helical enstrophy (sometimes called super- 
hclicity |40l), a pseudo-scalar as helicity itself. Note that, 
locally (as opposed to globally), helicity can be produced 
through alignment of vorticity and shear |4l| . However, 
globally and for an initially helical flow, Zh is respon- 
sible for the viscous decay of helicity. In the absence of 
dissipation, helicity is conserved for non stratified flows, 
while the first term on the r.h.s. of Eq. ((T^ can act as 
a source or a sink for stratified flows, thus breaking the 
conservation. 

In Fig. [1] (a) and (b) we first show the temporal evo- 
lution of the kinetic enstrophy, Zv{t), and potential en- 
strophy, Zp{t). All the flows have Fr sa 0.022, with, as 
initial conditions, either an ABC flow, an ABC2C flow, 
or a Taylor-Green flow. For comparison, one case with 
ABC initial conditions, is unstratified. The time in all 
the figures is expressed in units of the initial eddy turn- 
over time tnl = 27r/(fcoJ7o)- We see that the produc- 
tion of enstrophy, and therefore, the transfer of energy 
to small scales, is damped substantially in the presence 
of stratification, both for the ABC and TG initial condi- 
tions, although its maximum is not considerably delayed. 
This is expected because the effect of waves through the 
buoyancy forces is to reduce the nonlinear interactions, 
as well as to suppress, in part, the vertical velocity com- 
ponent. Observe that the stratified TG flow displays a 



significantly stronger peak for the kinetic and potential 
enstrophy than the stratified ABC-like flows for which 
the nonlinear terms are (initially) equal to zero. 

In Fig.[l](c) and (d) the temporal evolution of the total 
helicity and relative helicity for the same flows is shown. 
In the absence of stratification and for ABC initial condi- 
tions, helicity decays rather rapidly (exponentially after 
an initial nonlinear phase) and is close to zero for t > 6. 
On the contrary, stratified ABC and ABC2C flows dis- 
play slow decays of helicity, with strong oscillations at 
first when Vz{t = 0) 7^ 0, linked to gravity waves. Finally, 
we observe that the TG flow has zero initial helicity, and 
none is created by stratification. It is interesting that in 
the stratified cases helicity decays almost linearly with 
time, and much more slowly than the energy, which, as 
will be shown later and as is often the case for turbulent 
flows, decays as a power law with time. 

In the absence of stratification the relative helicity 
ay (t) decays rapidly, the flow becoming closer to mirror- 
symmetric. On the other hand, when gravity is switched 
on and initial helicity is non-zero, av(t) approaches the 
maximum value. This may be simply due to the fact 
that helicity decays more slowly than energy, leading to 
a growth of their ratio. More precisely, it can be seen 
from the data shown in Fig. [1] and later in Fig. [S] 
that dtH « {dtEvY''^{dtZvY''^, hence, in Eq. @ the 
denominator increases much faster than the numerator, 
causing cry to grow. 

The flows analyzed until now are well-ordered, cen- 
tered at large scales and with phase coherence between 
modes att = Q. We also examined initial conditions with 
randomized phases yet maintaining a high initial relative 
helicity (see Table HI)) . In Fig. [5] (a) and (b) the behavior 
of Zv(t) and Zp{t) for random flows with /cq = 3 — 4 or 
/co = 2 is shown. It is seen that the unstratified flow dis- 
plays a higher peak of the kinetic enstrophy, as noticed 
before, followed by the less stratified case {Fr « 0.044); 
while, when Fr « 0.022 the case with ko = 2 has a lower 
peak than the case with larger kg. This is likely due to 
the fact that when feg = 2 the flow preserves some helic- 
ity, which inhibits the energy decay. A similar behavior 
is seen for the potential enstrophy: the peak is higher 
for the less stratified case, and is lower when helicity is 
better preserved. 

In random flows, the time evolution of the helicity and 
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(a) (c) 




FIG. 1. Temporal evolution of the kinetic enstrophy (a), potential enstrophy (b), total helicity (c) and relative helicity (d), 
for several initial conditions, namely an unstratified ABC (run 6, ABCA^'O, dotted line), and three flows with Fr ~ 0.022: an 
ABC (run 7, ABC, solid line), an ABC2C (run 3, ABC2C, dash-dotted line), and a TG flow (run f, TG, dashed line). The 
helicity of the TG flow, initially zero, remains negligible (dashed line), whereas for the two other stratified and helical flows, 
the helicity decay is remarkably slow compared to the non-stratified case (dotted line). Note the quasi-maximal helical state 
reached in the ABC2C case (dash-dotted), especially clear in the evolution of (Jv{t)- 



the relative helicity reveals a more complicated behav- 
ior than in ABC and TG flows. It is indeed observed 
that if the injection wave number is fco = 2 (or fcp = 1, 
not displayed) , then, consistently with what we reported 
in the ABC cases, the relative helicity grows, while, for 
smaller-scale initial conditions, the relative helicity de- 
cays to a value close to zero (see Fig. [2]). If however, 
the Fr number is increased, a small growth of the rela- 
tive helicity is observed also for the fcg =3 — 4 case. A 
similar run at twice the Reynolds number (using a grid 
with a resolution of 512^ points) reveals no difference in 
the helicity behavior. The reason why random flows at 
fco = 2 behave similarly to ABC flows is likely due to the 
fact that their initial conditions resemble an ABC flow, 
especially with regard to the velocity distribution, with 
rather well-identified large scale structures; in the case 
of the ABC flow corresponding to three Beltrami waves, 
also called the Roberts flow, it resembles periodic cylin- 
drical vortices (see [i^] for a visualization of the ABC flow 
at t = Q). For higher initial fco, the randomness (and the 
existence of more modes that can be excited in a given 
Fourier shell) destroys these well-organized structures. 

Different snapshots of the x- and ^-components of the 
velocity, are given, respectively, in Fig. [3] and Fig. |4] for 
the ABC2Re and RND2Rc flows. Note that the horizon- 
tal cylindrical vortices, present in the ABC flow at i = 0, 
are being sheared and develop into a zig-zag shape at the 



peak of dissipation, reminiscent of the zig-zag instabil- 
ity found in laboratory experiments j43| . However, this 
does not appear as being as well defined in the random 
runs. At later times the x-component of the velocity is 
redistributed in horizontal layers of alternating sign; this 
structure is much better ordered for the ABC flow but 
is still visible for random flows. It is clear from inspec- 
tion of Fig. m that the vertical velocity keeps its cylindri- 
cal structuring in the ABC2Re flow, whereas it becomes 
more turbulent in the case of the RND2Re flow. It 
is also observed that, at a later stage, the ABC vertical 
columnar structures that survive alternate sign in time, 
suggesting a coupling with internal gravity waves. 

Finally, the effect of changing the Reynolds and Froude 
number on the time evolution of helicity for an ABC flow 
is shown in Fig. [S] By increasing the value of these 
dimensionless numbers, the total helicity increases ac- 
cordingly. The same phenomenon is observed, even if in 
smaller proportion, in the case of random flows. 



B. Cyclostrophic balance 

By measuring separately the flrst and second term on 
the r.h.s. of the dH/dt expression we flnd that, at early 
times, characterized by high amplitude gravity waves, the 
flrst term dominates. After the peak of enstrophy, the 
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FIG. 2. Temporal evolution of the kinetic enstrophy (a), potential enstrophy (b), total helicity (c) and relative helicity (d), for 
random initial conditions: an initially unstratified flow (run 12, RNDA'^0, dotted line), a flow with Fr ~ 0.022 (run 13, RND, 
dashed hue), a flow with ito = 2 and Fr f» 0.022 (run 17, RNDfc2, solid line), and a flow with Re = 6000 and Fr « 0.044 (run 
15, RND2Fr, dash-dotted). Note that all flows are close to maximum helicity initially (see Table The insets in (c) and (d) 
give a zoom on the early-time evolution of helicity. 



two terms are small and balanced and cause dH/dt k, 0. 
For the ABC2C flow with w^(i = 0) = 0, there is no 
need for radiation of gravity waves to adjust to large-scale 
balance and the two terms are small and comparable at 
all times. This may explain the large values of residual 
helicity at late stages. 

We observe that both the ABC and the random ko = 2 
flows initially excite smaller amplitude oscillations of the 
first term than the ko = 3 — 4 random flows do. These 
oscillations cease roughly at the peak of the enstrophy, 
when nonlinear coupling of waves takes over the linear 
phase. If the oscillation amplitude is large, H{t) goes to 
zero even before the peak of enstrophy (this is observed 
for the fco = 3 — 4 random flows), an irremediable situa- 
tion in the case of decaying flows, while, if the oscillation 
amplitude is small, H(t) will have a non-zero residual 
value at the peak of enstrophy, and from there on it will 
decay very slowly (ABC and fco = 2 random cases). 

Increasing the Reynolds number does not seem to mod- 
ify the behavior of H{t) before the peak of the enstro- 
phy, because the dynamical evolution in this first stage 
is plausibly wave-dominated. If the Froude number is 
varied instead, from Fr = 0.022 to Fr = 0.044, the ini- 
tial oscillations are reduced and now the random case 
with ko = 3 — 4 also sees a growth of the relative helic- 



ity (see Fig. [5]). We have to conclude that the growth 
of relative helicity, or similarly the slow decay of helicity, 
is non-monotonic with Froude number since in the case 
= 0, a decay of av{t) is observed. Another distinc- 
tion which may help interpreting the different behavior 
of helicity in the ABC flow and the random flow, is that 
while in ABC flows helicity is organized in large scale co- 
herent structures, in the structureless random flow it is 
not (except for fcp = 2). 

Wc are thus led to conjecture that, similarly to the 
geostrophic balance in rotating stratifled turbulence, an- 
other large-scale balance, linked to dissipation, can take 
place in helical flows in the presence of stratification 
alone. Such viscous cyclostrophic balance would be con- 
sistent with dH/dt « 0, as observed for the Froude num- 
bers considered. By neglecting the time derivative and 
the nonlinear term in Eq. (1), taking the curl of the re- 
duced equation and dotting it with velocity, we arrive at 
the following balance equation: 

A^u • V X (6* e^) = -vu ■ V'^uj = z/u • V x V x w. (13) 

Note that at this stage there is no space average being 
computed; however, upon space integration, this gives 
the relationship between vertical gradients of the temper- 
ature and the super helicity, as derived before, a balance 
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FIG. 3. (Color online) Distribution of the x (horizontal) velocity Vx for the ABC2i?e flow (left) and the RND2_Re flow (right) 
at times (from top to bottom) t — 0, t = t* (i.e., at the peak of dissipation), and t — 20. Light corresponds to low negative 
values, dark to high positive values. The vertical direction is indicated by the upward (blue) arrow. Observe the strong shear 
layers that have developed by the peak in enstrophy, and which later on begin to smooth out. 



FIG. 4. (Color online) Distribution of the vertical velocity for the ABC2_Re flow (left) and the RND2_Re flow (right) at times 
(from top to bottom) t = Q, t = t* , and t = 20. Light corresponds to low negative values, dark to high positive values. Observe 
the more turbulent flow in the random case. 
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FIG. 5. Temporal evolution of the total helicity for ABC flows 
with Re = 6000 and Fr = 0.022 (run 8, ABC2_Re, solid line), 
Re = 6000 and Fr = 0.044 (run 9, ABC2Fr, dashed-dotted 
line), and Re = 3000 and Fr — 0.022 for reference (run 7, 
ABC, dots). 



that may hold at large scales where the waves dominate 
the dynamics. This balance cancels the two terms on the 
r.h.s. of Eq. ([12]) thus explaining the slow decay of he- 
licity when it is initially at large scales. What happens 
with this balance when the buoyancy Reynolds number is 
increased will require runs at higher resolution (or lower 
stratification) than what is considered in the present pa- 
per, and will be the topic of a subsequent study (see Q 
for a large run where the Ozmidov scale is clearly re- 
solved). However, it is interesting to note that viscous 
forces arc known to play other important roles in strat- 
ified turbulence, such as when gravity waves break and 
alter the conservation of total energy and of potential 
vorticity (PV), leading to the creation of PV dipoles (see 



C. Energy decay 

We now turn our attention to the temporal evolution of 
energy. We display in Figs. [5] and [7] the decay of the total 
energy for a variety of initial conditions and Fronde num- 
bers (see captions for details). Computations are per- 
formed for roughly 30 turnover times (and thus, for 3000 
Brunt- Vai'ssala periods, for Fr = 0.01); some power-laws 
are added to guide the eye. When examining separately 
the temporal decay of the potential and kinetic energy, 
they evolve in similar ways but with strong oscillatory 
energy exchanges, while the oscillations due to gravity 
waves disappear when considering the total energy. 

In Fig. [6l it is striking to notice that the decay of en- 
ergy can be very different for different flows with the same 
external parameters (i.e., Reynolds and Fronde num- 
bers). If t* is the time at which dissipation sets in, that 
is, the maximum of enstrophy, the decay, in the absence 
of stratification, would follow a. ^ {t* — t)~^ law (dotted 
line) given that the growth of the integral scale is pre- 
vented by being in a "box-limited" case. Considering the 




FIG. 6. Temporal evolution of the total energy, with different 
scaling laws given as indications. ABC flow (run 7, ABC, solid 
line), ABC2C flow (run 3, ABC2C, dash-dotted line), and TG 
flow (run 1, TG, dashed line), the latter two having initially 
Vz — 0. All runs have Re ~ 3000 and Fr ~ 0.02 except for 
the ABC flow with A'^ = which is shown for comparison 
with a dotted line (run 6, ABCA'^0). Note the particularly 
slow decay for ABC runs 3 and 7. 

non-helical TG flow, the decay seems to follow a power- 
law ^ {t* — t)~^ , after an initial ideal (inviscid) phase. 
This result is expected on the basis of slowing-down of 
nonlinear interactions (leading to energy transfer and its 
decay) because of waves, and a similar power-law decay 
has already been observed for stratified flows [1^ , for ro- 
tating flows [i^, as well as for flows in the presence of 
a magnetic field [i^. For the ABC stratified fiow the 
decay of energy is substantially slower, with a power law 
of —1/3 as opposed to —1. On the other hand, we recall 
that helicity decays linearly and at a very slow rate, as 
shown in the preceding section. 

The slow decay of energy has been observed previously 
for rotating turbulence in the presence of helicity and the 
different power laws one may expect have been reviewed 
for a variety of cases takiiig into account the invariance of 
both energy and helicity [43|. A (t* — t)^^/"^ law is found 
on phcnomcnological grounds for helical rotating fiows 
based on the fact that helicity plays a role in the dynam- 
ics: it dominates the energy transfer to small scales and 
alters the spatial scaling laws for the energy spectra and 
higher order structure functions as well [4^, (see also 
[50| for review). The similar decay law found in this pa- 
per for ABC-like stratified fiows may be due to the fact 
that, as shown here, the helicity is quasi-conserved by 
the dynamics, due to a cyclostrophic balance, and thus 
leads to the same slow decay. We thus conclude that, for 
flows with either rotation or stratification, the presence 
of helicity considerably slows-down the energy decay, and 
measurably so, leading to persistent structures, whereas 
for the unstratified non-rotating case, helicity delays the 
onset of the decay but does not alter its rate of decrease 

In Fig.[7]wc examine the possible variation of the decay 
rate of energy with several factors: Froudc and Reynolds 
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FIG. 7. Temporal evolution of the total energy Erit), with different scaling laws given as indications for: (a) TG flows at 
different Froude numbers; Fr ^ 0.011 (run 2, TGFr/2, dashed line), and Fr ^ 0.022 (run 1, TG, dotted line); (b) ABC flows 
with Re = 6000 and Fr = 0.022 (run 8, ABC2i?e, dash-dotted line). Re = 6000 and Fr = 0.044 (run 9, ABC2i^r, solid line). 
Re = 3000 and Fr = 0.011 (run 11, ABCFr/2 dashed line), and Re = 3000 and Fr = 0.022 (run 7, ABC, dotted line); (c) 
random flows with Re = 6000 and Fr = 0.022 (run 14, RND2i?e, solid line). Re = 6000 and Fr = 0.044 (run 15, RND2fr, 
dash-dotted line), and Re — 3000 and Fr — 0.022 for reference (run 13, RND, dotted line). Note the significantly slower decay 
in (b) for the case of the ABC flow. 



numbers, as well as TG, ABC, and random initial con- 
ditions. We observe that for the TG flow, the decay is 
first faster (^ t~^) and then slows down; this could be 
due to the fact that, as time evolves and energy decays, 
the Froude number of a given fiow becomes substantially 
smaller; note that the later decay in Fig.[71[a) corresponds 
to an exponential regime, as can be easily seen display- 
ing the data in lin-log coordinates (not shown). If the 
Froude number is decreased the decay is delayed (dashed- 
line) , while still following the same power law. A similar 
effect is observed for ABC flows Fig. [TJb) where results 
for two different Reynolds numbers are reported: the 
higher the Reynolds number, the shallower the energy 
decay. 

For the random cases [Fig. El^c)], when /cq = 3 — 4, he- 
licity has become negligible by the peak of the enstrophy, 
as previously observed, and the decay follows an approx- 
imate law. When fcg = 2 there are too few modes 
to excite in building a random initial condition, and at 
early times the randomness of such flows is rather weak, 
with well-ordered columns in the initial conditions, rem- 
iniscent of the ABC flow itself (see, e.g.. Fig. [3]). Helicity 
then does not go to zero, and the decay follows, as in the 
ABC case, a —1/3 power law (this case is not shown). 

In general terms, it seems reasonable to conclude, for 
stratified flows, that if helicity vanishes by the time of 



occurrence of the peak of the enstrophy, the decay will 
follow a -1 power law, while if the initial helicity value is 
high and it survives up until the maximum of dissipation, 
the decay will follow a -1/3 power law. 

D. Energy spectra 

The effect of strong and large scale helicity on the tem- 
poral decay rate of the total energy is thus striking. We 
now address the related question of the possible role of 
helicity on energy distribution among Fourier modes. 

In Fig. [5] we show the isotropic (fc), parallel (fc||) and 
perpendicular {k±) total energy spectra, Et, together 
with the parallel and perpendicular potential energy 
spectra, Ep, for four ABC flows with different Froude 
and Reynolds numbers (see caption). The spectra are 
plotted at the peak of enstrophy, and are averaged over 
~ 0.3 — 0.4 eddy turnover times. As general remarks, 
we note that the energy is rather strongly peaked at fco, 
and we observe that, at high wavenumbers the spectra 
are more developed when the Froude number is higher, 
since nonlinear interactions are comparatively more effi- 
cient. The survival of the initially excited large scales is 
expected, since both wave interactions and helicity can 
inhibit the energy transfer to smaller scales. This fact 
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FIG. 8. Isotropic (a), parallel (b), and perpendicular (c) total energy spectra, and parallel (d) and perpendicular (e) potential 
energy spectra for ABC flows with Re = 6000 and Fr = 0.022 (run 8, ABC2i?e, solid line), Re = 6000 and Fr = 0.044 (run 9, 
ABC2Fr, dash-dotted line). Re = 6000 and Fr = 0.088 (run 10, ABC4Fr dashed line), and Re = 3000 and Fr = 0.022 (run 
7, ABC, dotted line). Slopes are indicated as a reference. Note that the rather flat spectrum in the total energy stems from 
vertical variations (left, b and d), whereas the perpendicular spectra (right, c and e) are close to a k'^^ scaling for all values of 
parameters. 



may explain the rather shallow spectrum at intermedi- 
ate wavenumbers, followed by the very steep spectrum 
at high wavenumbers, observed in Fig. [5] (a) and (b). 
In the isotropic energy spectrum [Fig [5] (a)], the transi- 
tion between the rather flat behavior and the steep high- 
wavenumber scaling law is marked by a 'knee' around 



sa 30, which is possibly coincident with the buoyancy 
wavenumber ki,: ki, w 25 for Fr — 0.022, and ki, w 13 
for Fr = 0.044. The same knee is observed in the total 
and potential parallel energy spectra, being so entirely 
attributable to the behavior of the flow with respect to 
variations in the vertical direction, a fact that strength- 
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FIG. 9. Isotropic (a), parallel (b), and perpendicular (c) helicity spectra for the same flows as in Fig. [8l run 8 (solid line), run 
9 (dash-dotted line), run 10 (dashed line) and run 7 (dotted line). Note the flat spectra at large scales, stemming again from 
vertical variations. 



ens its identification with kf,. 

The perpendicular spectra decay consistently with a 
power law, as found in several numerical studies for 
this parameter regime 0, [5l| . 

Flat spectra have already been documented in the lit- 
erature for stratified turbulence and various interpreta- 
tions have been provided. A flat energy spectrum has 
been derived in [2^1 on the basis of a superposition of 
independent horizontal layers, each with a Gaussian pro- 
file, as a model of the vertically-sheared horizontal lay- 
ers associated with the zig-zag instability [1^ (as well as 
with geostrophic balance when one takes rotation into ac- 
count). Indeed, the thickness of these layers is thought to 
be the buoyancy scale, as advocated in [i^ using a scal- 
ing argument that states that the Froude number based 
on a vertical length scale has to be of order unity. This 
has already been studied in a series of direct numerical 
simulations (see, e.g., [13, [13 for more details). Here 
we hypothesize that flat spectra may be related to the 
presence of high helicity values as we discuss below com- 
menting on Fig. [S] 

The knee we have identified in the evolution of energy 
spectra as a function of fc|| is also not a new feature of 
stratified turbulence. It has been attributed to a change 



of regime for the total energy (or its horizontal compo- 
nent) at high resolution (specifically with the dissipation 
wave number sufficiently resolved) at the buoyancy scale 
[53| . On the other hand, in [5l|, both the wave mode 
(corresponding to the vertical velocity) and the vortex 
mode (corresponding to the horizontal velocity) display 
a knee in terms of fc|| for highly stratified flows, identified 
as the buoyancy scale. Our data leads us to argue (see 
below) that another, rather independent phenomenon, 
may be happening as well, namely the role played by the 
helicity distribution among modes. 

The isotropic, parallel and perpendicular helicity spec- 
tra are shown in Fig. |9]for the same ABC flows (note that 
the absolute value of helicity is displayed) . Noteworthy, 
similarly to the total and parallel isotropic energy spectra 
are the following points: (i) the spectrum is fiat at large 
scale; and (ii) a knee is present for the total and parallel 
helicity; in addition we also observe a secondary change 
in the spectra at a smaller scale which seems to coincide 
with a change of sign of helicity. Since, according to the 
Schwartz inequality we have H{k) ^ kE{k), and since 
helicity is large for these flows, it is reasonable to think 
that the behavior of the energy spectra is influenced by 
the behavior of helicity. 
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FIG. 10. Isotropic (a), parallel (b), and perpendicular (c) total energy spectra, and parallel (d) and perpendicular (e) potential 
energy spectra for random flows with Re = 6000 and Fr = 0.022 (run 14, RND27?e, solid line). Re = 6000 and Fr = 0.044 (run 
15, RND2Fr, dash-dotted line). Re = 6000 and Fr = 0.088 (run 16, RND4fr, dashed line), and Re = 3000 and Fr = 0.022 
(run 13, RND, dotted line). Slopes are indicated as a reference. No knee in the spectra is discernible, contrary to ABC flows. 



As the temporal data indicates that the ABC flow be- 
haves differently from random flows, we now examine 
the spectral behavior for the random case. In Fig. [TU] 
the energy spectra for random flows centered initially on 
wavenumber ko — 3, 4 are presented. We recall that at 
the peak of enstrophy, the residual helicity is close to zero 
for these runs. In these cases there is no clear signature 
of a knee or an abrupt slope change for the energy. The 



higher the Re and Fr numbers, the shallower the spectra 
at small scales. The peak at fco is also noticeably reduced 
with respect to the ABC cases, indicating more efficient 
energy transfer from the energy-containing scales. 

Finally, in Fig. [11] the helicity spectra for the random 
cases initially centered at fco = 3, 4 are shown. The value 
of helicity is much smaller than for the ABC runs, as 
noted in Table HIl the residual relative helicity ay at the 
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FIG. 11. Isotropic (a), parallel (b), and perpendicular (c) helicity spectra for the same random flows as in Fig. [10] run 14 
(solid line), run 15 (dash-dotted line), run 16 (dashed line) and run 13 (dotted line). 



peak of enstrophy is close to 2%, whereas for the ABC 
runs it is close to 20%. At the end of the runs, these 
values are respectively ~ 11% and « 80%, i.e., almost as 
strong as it can be in the latter case. A slope change and 
a knee arc still identifiable at a wavcnumbcr again com- 
parable to the buoyancy wavenumber fcb, but the helicity 
is now too weak to influence the behavior of the energy 
spectrum. Note also that at scales larger than the buoy- 
ancy scale (but smaller than the initial condition scale), 
the helicity spectrum is again flat, likely due to the lay- 
ered structure in the vertical, whereas the perpendicular 
spectrum is steep. 

To conclude, it seems likely, for helical flows here rep- 
resented only by ABC initial conditions, that it is the 
helicity which is responsible for the change of slope in the 
energy spectrum close to the buoyancy scale, and not the 
other way round. Assuming that the Schwarz inequality 
is fulfilled, a flat helicity spectrum leads to E{k) ^ k~^, 
compatible with what we observe for the energy in the 
ABC case. A similar result obtains for the random ini- 
tial condition centered on fcg =2. This point will need 
further studies at high resolution so that different energy 
and helicity balance can take place in flows with sufficient 
scale separation. 



IV. CONCLUSION 

We have performed numerical simulations of stratified 
turbulence and have shown that initial conditions with 
sufficiently strong large-scale helicity produce measur- 
able effects. This is observable in both the energy de- 
cay, which is found to be substantially slower than in the 
non-helical case, and in the energy (and helicity) spec- 
tra, which display a flat distribution with a marked knee 
in the vicinity of the buoyancy scale Lh] beyond L;,, the 
spectra are signiflcantly steeper for the cases examined 
in this paper, which are strongly stratifled, with buoy- 
ancy Reynolds number of order unity and up to 46 (see 
Table |l| . The knees are clearly identifled with variations 
of these spectra with the vertical wavenumber, and the 
knees seem to be close to the buoyancy scale, in agree- 
ment with the interpretation of Lh as the scale at which 
Fr - 1. 

Thus, for helical stratifled turbulence, the energy de- 
cay is slowed down by the presence of helicity in a way 
reminiscent to that observed in the rotating case. How- 
ever, it should be noted that helicity is an ideal invariant 
of rotating flows, whereas it is not a priori conserved for 
ideal stratifled flows. While in the former case the slow 
down in the decay is associated with a direct cascade 
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of the helicity, in the latter case a very slow (linear in 
time) decay of helicity (resulting in a quasi-conservation) 
can be attributed to a cyclostrophic balance where dis- 
sipation and gravity equilibrate at large scales once the 
nonlinear terms start transferring energy to small scales 
after an initial wave-dominated regime. These results 
are obtained in cases where the evolution of large scales 
is somewhat inhibited, because the initial scale of the 
flow is comparable to the overall size of the computa- 
tional box, the so-called box-limited conditions in which 
the growth of the integral scale is blocked. 

Moreover, in highly turbulent flows such as those found 
in geophysics and astrophysics, the eddy turn-over time 
becomes shorter than tw at smaller scales, for a suffi- 
ciently large Reynolds number. In that case, isotropic 
turbulence may recover, beyond the so-called Ozmidov 
scale Loz ^ [e/N^Y^^. Assuming further that this scale 
is larger than the dissipation scale idiss ~ [e/i^^]^^^ leads 
to the condition that the buoyancy Reynolds number 
TZ = ReFr^ be larger than unity. It was found recently 
that power laws of the energy spectra scale with TZ, as 
opposed to the Froude number itself [s^. The condi- 
tion 7?, >> 1 is easily fulfilled in geophysical flows but 
difficult to realize in numerical simulations, where the 
highest to date is found in realized on a grid of 



4096^ * 1024 points; similarly, laboratory experiments 
function at moderate TZ [s^. Further studies of that 
regime are necessary but will require substantial numer- 
ical resources. 
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